Figure 4 and supplementary figure 4

Author

Niccolò Arecco

1 Intro

We performed histone PTMs DIA MS and I analysed the raw files with EpiProfile v2.1 (Yuan et al. 2018). Here I parse the output and analyse the PTMs values with DEP.

Acknowledgement: I’d like Zuo-Fei Yuan and Simone Sidoli from Ben Garcia lab for their help with data analysis and sample preparation respectively.

2 Set Up

To fetch files on the CRG cluster I mount the server on my local computer with the following command in the terminal.

Code
sshfs narecco@ant-login.linux.crg.es:/users/mirimia ~/mnt

2.1 Packages

Load packages required for the analysis and suppress any message. Check the Section 5 section at the end for package versions.

Code
library(dplyr, warn.conflicts = F, quietly = T)
library(readr)
library(readxl)
library(DESeq2)
library(csaw)
library(pheatmap)
library(DEP)
library(Cairo)
# library(knitr)
library(DT)

Load my R package niar to use my custom made function to parse vast-tools output.

Code
if ( ('niar' %in% .packages(all.available = TRUE)) == TRUE ) {
  library('niar')
} else if ( ('niar' %in% .packages(all.available = TRUE)) == FALSE ){
  message("The package niar is not installed so I'm trying to load it from ",
          "the local repository of the package")
  if ( dir.exists('~/mnt/narecco/software/R/niar' ) )  {
    devtools::load_all(path = '~/mnt/narecco/software/R/niar')
  } else {
    stop("Can't find the local repo of the niar package! ",
         "You must install it with:\n",
         "devtools::install_github('Ni-Ar/niar') ")
  }
} else{
  stop("Can't understand if 'niar' package was installed beforehand")
}
The package niar is not installed so I'm trying to load it from the local repository of the package
ℹ Loading niar

2.2 Functions

Define a plot style.

Code
theme_classic(base_size = 6.3, base_family = "Arial") +
  theme(plot.title = element_text(size = 8, vjust = -1, hjust = 0),
        plot.background = element_blank(),
        strip.background = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(linewidth = 0.1),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 5),
        axis.text = element_text(colour = "black"),
        axis.text.x = element_text(margin = margin(t = 0, unit = "mm")),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_line(linewidth = 0.1),
        axis.ticks.length.y = unit(1, "mm"),
        legend.margin = margin(l = 0, t = 0, unit = "mm"),
        legend.box.margin = margin(l = 0, unit = "mm"),
        legend.box.background = element_blank(),
        legend.box.spacing = unit(0.5, "mm"),
        # legend.title = element_blank(),
        legend.text = element_text(margin = margin(l = -10, unit = "mm")),
        legend.text.align = 0,
        legend.key.height = unit(2, "mm"),
        legend.key.width = unit(7.5, "mm"),
        legend.key = element_blank(),
        legend.background = element_blank(),
        legend.spacing.x = unit(3, 'mm')
        ) -> alluvial_thm 

I modify the geom_stratum() function just to add a very small linewidth the the black boders of the colums. With linewidth = .05, I don’t have to change them manually later

Small function to parse the PTM abundance as 10 to the power scientific notion.

Code
scientific_10 <- function(x) {
  parse(text=gsub("e\\+", " %*%10^", scales::scientific_format()(x)))
}

2.3 Directories & File Paths

Here I organise all the file and directories paths I need to run the analysis and define where to save the processed tables and figures.

Code
oneDrive_Dir <- file.path("~/OneDrive - CRG - Centre de Regulacio Genomica/Suz12_AS_project")  
code_dir_fig <- file.path(oneDrive_Dir, "_Code/Fig4")
tbl_dir_fig4 <- file.path(code_dir_fig, "tables")
pdf_dir_fig4 <- file.path(code_dir_fig, "pdfs")

if (!dir.exists(pdf_dir_fig4)) { dir.create(pdf_dir_fig4, recursive = T) }
if (!dir.exists(tbl_dir_fig4)) { dir.create(tbl_dir_fig4, recursive = T) }
Code
metadata_path <- file.path(tbl_dir_fig4, 'Histone_PTMs_metadata.xlsx')

WT_Dex4_histone_ratio_path <- list.files(tbl_dir_fig4,
                                         pattern = "FullSet1_Run2_DIA_histone_ratios_Suz12dExon4.xls",
                                         full.names = T)

Resce_KO_histone_ratio_path <- list.files(tbl_dir_fig4,
                                          pattern = "FullSet2_Run1_DIA_histone_ratios_Suz12KO.xls",
                                          full.names = T)

stopifnot(file.exists(metadata_path))
stopifnot(file.exists(WT_Dex4_histone_ratio_path))
stopifnot(file.exists(Resce_KO_histone_ratio_path))

3 Main Figure Panels

Import metadata.

Code
mtdt <- read_excel(metadata_path) |> 
  setNames(c("Sample", "Short_Name", "Clone_Name", "Condition", "Replicate"))

mtdt$Condition <- factor(mtdt$Condition, 
                         levels = c("WT", "CSexon4", "Dexon4", "KO", "KO_L", "KO_S") )

3.1 H3 K27-K40 peptide modification ratios

Import DIA data quantification using automatic EpiProfile v2.1 quantifications.

Code
read_EpiProfile_histone_ratio(WT_Dex4_histone_ratio_path) |>
  left_join(y = mtdt, by = "Sample") |>
  mutate(Facility_ID = Sample) |>
  mutate(Sample = Short_Name) |>
  tidy_hPTMs(split_double_PTMs = T) -> dia_auto_set1

read_EpiProfile_histone_ratio(Resce_KO_histone_ratio_path) |>
  left_join(y = mtdt, by = "Sample") |>
  mutate(Facility_ID = Sample) |>
  mutate(Sample = Short_Name) |>
  tidy_hPTMs(split_double_PTMs = T) -> dia_auto_set2

data_dia_auto <- rbind(dia_auto_set1, dia_auto_set2)

Tidy up the information for histone peptide H3 K27- K36. Here I split the H3 K27 and K36 signal coming form the same peptide. In this way I can pile up the all methyl-forms.

Code
data_dia_auto |>
  subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
  mutate(K27 = case_when( is.na(Sub_PTM1) & grepl(pattern = "^K27", x = Modification) ~ Modification,
                          grepl(pattern = "unmod", x = Modification) ~ "unmod" ) ) |>
  mutate(K36 = case_when( is.na(Sub_PTM2) & grepl(pattern = "^K36", x = Modification) ~ Modification, 
                          grepl(pattern = "unmod", x = Modification) ~ "unmod" ) ) |>
  relocate(K27, K36, .after = Sub_PTM2) |>
  unite(col = "K27", c("Sub_PTM1", "K27"), sep = "", na.rm = T, remove = T ) |>
  unite(col = "K36", c("Sub_PTM2", "K36"), sep = "", na.rm = T, remove = T ) |>
  mutate(K27 = case_when(grepl(pattern = "^$", x = K27) ~ "unmod",
                         grepl(pattern = "[aA-zZ]", x = K27) ~ K27 )) |>
  mutate(K36 = case_when(grepl(pattern = "^$", x = K36) ~ "unmod",
                         grepl(pattern = "[aA-zZ]", x = K36) ~ K36 )) |>
  # When one mark is not found I label it unmodified ('unmod') but it could also
  # be labelled as "only the other mark".
  # mutate(K27 = case_when(grepl(pattern = "^$", x = K27) ~ "Only\nK36me",
  #                        grepl(pattern = "[aA-zZ]", x = K27) ~ K27 )) |>
  # mutate(K36 = case_when(grepl(pattern = "^$", x = K36) ~ "Only\nK27ac-me",
  #                        grepl(pattern = "[aA-zZ]", x = K36) ~ K36 )) |>
  relocate(K27, K36, .after = PTM) |>
  mutate(Histone = factor(Histone, levels = c("H3", "H3.3") ) ) -> data_K27_36

I can also deconvolve the K27 and K36 signals from the same peptide and pile-up each signal.

Code
data_K27_36 |>
  group_by(Condition, K27, Histone) |>
  mutate(Average_Ratio = mean(Ratio)) |>
  mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) )|>
  select(Histone, Condition, Average_Ratio, K27, Short_Name, Ratio) |>
  mutate(Ratio_Sum = sum(Average_Ratio) / 3) |>
  mutate(Percent = round(Average_Ratio*100, 2) ) |>
  select(Condition, Histone, K27, Average_Ratio, Ratio_Sum, Percent) |>
  unique() |>
  mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
  arrange(Condition, Histone, K27) |> ungroup() -> data_K27

Same as above but for K36.

Code
data_K27_36 |>
  group_by(Condition, K36, Histone) |>
  mutate(Average_Ratio = mean(Ratio)) |>
  mutate(K36 = factor(K36, levels = c("K36me3", "K36me2", "K36me1", "unmod") ) ) |> 
  select(Histone, PTM, Condition, Average_Ratio, K36, Short_Name, Ratio) |>
  mutate(Ratio_Sum = sum(Average_Ratio) / 3) |>
  mutate(Percent = round(Average_Ratio*100, 2) ) |>
  select(Condition, Histone, K36, Average_Ratio, Ratio_Sum, Percent) |>
  unique() |>
  mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
  arrange(Condition, Histone, K36) |> ungroup() -> data_K36

For specific histone PTMs one can explore the results in the following tables.

Overview of H3 and H3.3 K27 acetylation and methylation quantifications

Code
datatable(data_K27, rownames = FALSE, filter = 'top', 
          options = list(pageLength = 5, autoWidth = TRUE) )

Overview of H3 and H3.3 K36 methylation quantifications

Code
datatable(data_K36, rownames = FALSE, filter = 'top', 
          options = list(pageLength = 5, autoWidth = TRUE) )

3.1.1 H3 K27

Plot figure 4B panel with the pileup of figure H3K27 methylation and acetylation.

Code
data_K27_36 |>
  group_by(Condition, K27, Histone) |>
  mutate(Average_Ratio = mean(Ratio)) |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3) |> # 3 replicates
  subset(Histone %in% "H3") |>
  mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) ) |>
  ggplot() +
  aes(x = Condition, stratum = K27, alluvium = K27, 
      y = Ratio_Sum, fill = K27) +
  geom_alluvium(alpha = 0.5, width = 0.4) +
  geom_stratum2(alpha = 1, width = 0.6, colour = "black") +
  geom_fit_text(aes(label = K27 ),
                stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"),
                padding.y = grid::unit(0.2, "mm"),
                size = 7, family = "Arial", show.legend = F) +
  scale_fill_manual(values = c("K27ac" = "#feac81", "K27me1" ="#ced1af", 
                               "K27me2" = "#748f46", "K27me3" = "#47632a", 
                               "unmod" = "gray" ),
                    name = 'H3') +
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 
                              'KO+L', 'KO+S')) +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing 
  labs(y = "Average Ratio") +
  alluvial_thm -> p_Alluvial_K27

p_Alluvial_K27 

Save to pdf.

Code
ggsave(filename = "Fig4B_H3_only_K27_Mrgd_Alluvial.pdf", path = pdf_dir_fig4, 
       plot = p_Alluvial_K27, device = cairo_pdf, units = "cm",
       width = 6.5, height = 4.85)

3.2 Differential abbundance analysis of histone PTMs with DEP

Use the area under the peak as counts to be analysed using DEP framework.

Code
rbind( read_EpiProfile_histone_ratio(WT_Dex4_histone_ratio_path),
       read_EpiProfile_histone_ratio(Resce_KO_histone_ratio_path) ) |>
  left_join(y = mtdt, by = "Sample") |>
  mutate(Facility_ID = Sample) |>
  mutate(Sample = Short_Name) -> data_dia_auto_not_tidy

One small caveat is that the areas under the peaks are very large numbers or at least larger than 2147483647 which is the highest number R can store as integer. So to convert the areas to integers from numeric type one could use the package bit64 that stores integers in 64 bits instead of the default 32 bit. My solution was however simpler, I divided the all table by 1000 to scale down by 3 order of magnitude. I then convert the dataframe to matrix.

I tried to add a pseudocounts to the data but it actually doesn’t help with the analysis and is better to input the areas as they are. The function make_unique is the first step of the analysis of DEP and creates a unique data frame.

Code
data_dia_auto_not_tidy |>
  tidy_hPTMs() |>
  select(First_Col, PTM, Sample, Area) |>
  # reduce all areas by a 1000 to avoid having to coerc to integer big numbers
  mutate(Area = Area / 1e4) |>
  pivot_wider(id_cols = c(First_Col, PTM), names_from = Sample, values_from = Area) |>
  # Coerce the areas to integers
  mutate(across(.cols = ends_with(c("_1", "_2", "_3")), .fns = as.integer )) |>
  make_unique(names = "First_Col", ids = "PTM") -> data_unique

Create a metadata compatible for DEP set the KO as the reference sample as I’m interesting in the rescue efficiency relative to the KO.

Code
# Set metadata colnames to be as `make_se` wants them to be 
colnames(mtdt) <- c("Sample", "label", "Clone_Name", "condition", "replicate")

mutate(mtdt, condition = factor(condition, levels = c("WT", "CSexon4", "Dexon4",
                                                      "KO", "KO_L", "KO_S"))) |>
  print(n = 1 ) -> mtdt_dep
# A tibble: 18 × 5
  Sample                          label Clone_Name condition replicate
  <chr>                           <chr> <chr>      <fct>         <dbl>
1 2022MZ006_IVMO_001_02_25pto_DIA WT_1  ZWT        WT                1
# ℹ 17 more rows

Find the columns that contain the samples areas in the data_unique dataframe and create a SummarizedExperiment object to be used in for the DEP analysis. The first step is to do a log2 transformation of the area counts.

Code
Samples_area_columns <- grep(pattern = paste(mtdt_dep$label, collapse = "|"), colnames(data_unique))
data_se <- make_se(data_unique, columns = Samples_area_columns, expdesign = mtdt_dep)

Now the data_se contains the modifications in the rows (assuming they are proteins) and in the column the area.

Code
data_filt <- filter_missval(data_se, thr = 0)

Check PTMs frequency across samples.

Code
plot_frequency(data_filt)

Overall most modifications are found in all samples.

Code
data_norm <- normalize_vsn(data_filt)

Verify the fit.

Code
meanSdPlot(x = data_norm)

Seems okay.

Check normalisation effect on data.

Code
plot_normalization(data_se, data_filt, data_norm) +
  scale_fill_manual(values = c("WT" = "royalblue3", "Dexon4" = "goldenrod",
                               "CSexon4" = "#B0461C", "KO" =  "gray45", 
                               "KO_L" = "mediumpurple3", "KO_S" = "darkorange1")) +
  theme(axis.text = element_text(color = 'black'),
        plot.background = element_blank(), 
        panel.background = element_blank())

Check missing values.

Code
plot_missval(data_norm)

I don’t see obvious missing values trends, besides only that sample 3 of CSexon 3 has less PTMs detected. I don’t think it’ necessary to impute missing values.

Code
# Test every sample versus control
data_diff <- test_diff(data_norm, type = "control", control = "WT")
Warning in test_diff(data_norm, type = "control", control = "WT"): Missing
values in 'data_norm'
Tested contrasts: CSexon4_vs_WT, Dexon4_vs_WT, KO_vs_WT, KO_L_vs_WT, KO_S_vs_WT
Warning: Partial NA coefficients for 5 probe(s)
Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There may be
too few input test statistics for reliable FDR calculations!

Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There may be
too few input test statistics for reliable FDR calculations!

Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There may be
too few input test statistics for reliable FDR calculations!

Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There may be
too few input test statistics for reliable FDR calculations!

Warning in fdrtool::fdrtool(res$t, plot = FALSE, verbose = FALSE): There may be
too few input test statistics for reliable FDR calculations!

Significance thresholds:

  • P-value <= 0.05

  • |Log2 Fold Change| >= 1.5

Code
dep <- add_rejections(data_diff, alpha = 0.05, lfc = 1.5 )

Check sample clustering a Gower’s distance matrix.s

Code
plot_dist(dep, significant = TRUE, pal = "PuOr") 

Check sample clustering.

Code
plot_heatmap(dep, type = "centered", kmeans = F, clustering_distance = "gower",
             col_limit = 4, show_row_names = TRUE,
             indicate = c("condition", "replicate"))
Warning: Missing values in 'dep'. Using clustering_distance = 'gower'

Get results.

Code
data_results <- get_results(dep)
data_results |>
  as_tibble() |>
  select(-significant) |>
  pivot_longer(cols = !c("name", "ID"),
               names_to = c("contrast", "Feature"),
               names_sep = c("_(vs_WT_|centered)"),
               values_to = "Value"
               ) |>
  mutate(Feature = ifelse(test = Feature == "", yes = "centered", no = Feature)) |>
  pivot_wider(names_from = "Feature", values_from = "Value") |>
  dplyr::select(-ID) -> res_all_tidy

colnames(res_all_tidy)[colnames(res_all_tidy) == "name"] <- 'First_Col'

Extract WT sample area of each PTM, to be used for the volcano plot.

Code
data_dia_auto |>
  subset(Condition == "WT") |>
  group_by(First_Col) |>
  mutate(WT_Mean_Area = mean(Area) + 1 ) |>
  relocate(WT_Mean_Area, .after = PTM) |>
  select(c(First_Col, WT_Mean_Area)) |>
  unique() -> WT_area

Get the PTMs info columns.

Code
histone_ptms_anno <- unique(data_dia_auto[, 1:10])

Add all the info the the tidy results dataframe.

Code
res_all <- left_join(x = res_all_tidy, y = histone_ptms_anno, by = join_by(First_Col))
res_all <- left_join(res_all, WT_area, by = join_by(First_Col) )

res_all |>
  mutate(Histone = factor(Histone, levels = sort(unique(res_all$Histone)) ) ) |>
  arrange(desc(Histone), Peptide_Start, Peptide_End) |>
  ungroup() -> res_all

Export all tidy results for supplementary table 3 as tab delimited, or as csv with UTF-8 Byte order mark which indicates to Excel the csv is UTF-8 encoded.

Code
write_delim(x = res_all, file = file.path(tbl_dir_fig4, 'DEP_analysis_histone_MS.tab'), 
            delim = '\t', col_names = T, append = F, quote = 'none', na = "NA",
            progress = F, escape = 'none')

write_excel_csv(x = res_all, 
                file = file.path(tbl_dir_fig4, 'DEP_analysis_histone_MS.csv'),
                 delim = ',', col_names = T, append = F, quote = 'none', 
                na = "NA", progress = F, escape = 'none')

Interactively explore the table

Code
res_all |>
  mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
  datatable(rownames = FALSE, filter = 'top', 
            options = list(pageLength = 3, autoWidth = TRUE) )

3.2.1 PTMs centred abundances

Plot centred abundance

Code
res_all |>
  subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
  subset(!is.na(centered)) |>
  ggplot(aes(x = centered, y = PTM, fill = contrast)) +
  facet_grid(rows = vars(Modification), scales = 'free' ) +
  geom_vline(xintercept = 0, linewidth = 0.4) + 
  geom_col(width = 0.75, 
           position = position_dodge(width = 0.75, preserve = 'single')) +
  geom_point(aes(size = WT_Mean_Area),
             shape = 21, position = position_dodge(width = 0.75),
             stroke = 0.2) +
  labs(x = "Centered PTM abundance") +
  scale_x_continuous(limits = c(-5, 5), oob = scales::oob_squish ) +
  scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
                  labels = scientific_10, name = "WT PTM\nabundance") +
  scale_fill_manual(values = c("WT" = "royalblue3", "Dexon4" = "goldenrod",
                               "CSexon4" = "#B0461C", "KO" =  "gray45", 
                               "KO_L" = "mediumpurple3", "KO_S" = "darkorange1"),
                    name = "Sample") +
  theme_classic() +
  theme(strip.text = element_blank(),
        strip.background = element_blank(),
        axis.text = element_text(colour = "black"),
        axis.title.y = element_blank(),
        legend.key.size = unit(1, 'mm'),
        panel.grid.major.y = element_line()) -> p_centered_K27_K36

p_centered_K27_K36

3.2.2 Volcano plot ∆ex4 vs WT

Prepare data for plot.

Code
res_Dex4 <- subset(res_all, contrast == "Dexon4" & !is.na(p.val))

res_Dex4 <- res_Dex4 |>
  mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
                              p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
                              p.adj > 0.05 ~ 'None') ) 

Plot figure 4C.

Code
ggplot(res_Dex4 ) +
  aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
  geom_point(shape = 21, stroke = 0.2) +
    annotate(geom = "label", x = 5.5,  y = 0.10, label = "∆ex4", 
           colour = "black", fill = "goldenrod", size = 2,
           label.padding = grid::unit(0.5, "mm"), 
           label.r = unit(0.25, "mm"), family = "Arial",
           label.size = grid::unit(0.125, "mm") ) +
  annotate(geom = "label", x = -5.5, y = 0.10, label = "WT", 
           colour = "white",  fill = "royalblue3", size = 2, 
           label.padding = grid::unit(0.5, "mm"), 
           label.r = unit(0.25, "mm"), family = "Arial",
           label.size = grid::unit(0.125, "mm") ) +
  # points on the right side
  geom_label_repel(data = subset(res_Dex4, p.adj < 0.09 & ratio > 0 ), 
                   aes(label = PTM), fill = 'white',
                   seed = 16, show.legend = F, segment.curvature = -1e-20, 
                   family = "Arial", size = 1.85, nudge_x = -0.25,
                   segment.color = 'black',verbose = F, 
                   box.padding = grid::unit(1, "mm"),
                   point.padding = grid::unit(0.55, "mm"),
                   label.padding = grid::unit(0.5, "mm"), 
                   label.size = 0.120, 
                   max.overlaps = 20)  +
  # points on the left side
  geom_label_repel(data = subset(res_Dex4,  p.adj < 0.1 & ratio < 0 ), 
                 aes(label = PTM), fill = 'white',
                 seed = 16, show.legend = F, segment.curvature = -1e-20, 
                 family = "Arial", size = 1.85, nudge_x = 0.25,
                 segment.color = 'black',verbose = F, 
                 box.padding = grid::unit(1, "mm"),
                 point.padding = grid::unit(0.55, "mm"),
                 label.padding = grid::unit(0.5, "mm"), 
                 label.size = 0.120, 
                 max.overlaps = 20)  +
  labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
  scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
  guides(size = guide_legend(override.aes = list(fill = "white"))) +
  scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-6.4, 6.4), n.breaks = 7) +
  scale_y_continuous(expand = expansion(add = c(0, 0.25), mult = 0), limits = c(0, NA)) +
  scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
                  labels = scientific_10, name = "WT PTM\nabundance") +
  theme_classic(base_size = 6, base_family = "Arial") +
  theme(panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_line(linewidth = 0.2),
        legend.position = c(0.91, 0.87), 
        legend.key = element_blank(),
        legend.key.size = unit(1.0, units = 'mm'),
        legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
        legend.direction = "vertical",
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.box.spacing = unit(-3, 'mm'),
        axis.text = element_text(colour = 'black'),
        axis.ticks.length = unit(1, units = 'mm'),
        axis.ticks = element_line(colour = 'black'),
        axis.title = element_text(size = 5),
        axis.line = element_line(linewidth = 0.2),
        plot.background = element_blank() ) -> p_Dex4

p_Dex4

Save to pdf.

Code
ggsave(path = pdf_dir_fig4, filename = "Fig4C_Dex4_vs_WT_Volcano_w_labels.pdf",
       plot = p_Dex4, device = cairo_pdf, units = "cm",
       height = 5.4, width = 7.0)

3.2.3 Other volcano plots

Prepare data of CSex4 vs WT

Code
res_CSex4 <- subset(res_all, contrast == "CSexon4" & !is.na(p.val))

res_CSex4 <- res_CSex4 |>
  mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
                              p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
                              p.adj > 0.05 ~ 'None') ) 

Note: nothing is differentially expressed in this contrast

Code
table(res_CSex4$Direction)

None 
 159 

Prepare data of KO vs WT

Code
res_KO <- subset(res_all, contrast == "KO" & !is.na(p.val))

res_KO <- res_KO |>
  mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
                              p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
                              p.adj > 0.05 ~ 'None') ) 

Plot KO vs WT

Code
ggplot(res_KO ) +
  aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
  geom_point(shape = 21, stroke = 0.2) +
    annotate(geom = "label", x = 8.5,  y = 0.10, label = "KO", 
           colour = "white", fill = "gray16", size = 2,
           label.padding = grid::unit(0.5, "mm"), 
           label.r = unit(0.25, "mm"), family = "Arial",
           label.size = grid::unit(0.125, "mm") ) +
  annotate(geom = "label", x = -12.5, y = 0.10, label = "WT", 
           colour = "white",  fill = "royalblue3", size = 2, 
           label.padding = grid::unit(0.5, "mm"), 
           label.r = unit(0.25, "mm"), family = "Arial",
           label.size = grid::unit(0.125, "mm") ) +
  # points on the right side
  geom_label_repel(data = subset(res_KO, p.adj <= 0.15 & ratio > 0 ), 
                   aes(label = PTM), fill = 'white',
                   seed = 16, show.legend = F, segment.curvature = -1e-20, 
                   family = "Arial", size = 1.85, nudge_x = -0.25,
                   segment.color = 'black',verbose = F, 
                   box.padding = grid::unit(1, "mm"),
                   point.padding = grid::unit(0.55, "mm"),
                   label.padding = grid::unit(0.5, "mm"), 
                   label.size = 0.120, 
                   max.overlaps = 20)  +
  # points on the left side
  geom_label_repel(data = subset(res_KO, p.adj <= 0.15 & ratio < 0 ), 
                 aes(label = PTM), fill = 'white',
                 seed = 16, show.legend = F, segment.curvature = -1e-20, 
                 family = "Arial", size = 1.85, nudge_x = 0.25,
                 segment.color = 'black',verbose = F, 
                 box.padding = grid::unit(1, "mm"),
                 point.padding = grid::unit(0.55, "mm"),
                 label.padding = grid::unit(0.5, "mm"), 
                 label.size = 0.120, 
                 max.overlaps = 20)  +
  labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
  scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
  guides(size = guide_legend(override.aes = list(fill = "white"))) +
  scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-15, 10), n.breaks = 7) +
  scale_y_continuous(expand = expansion(add = c(0, 0.25), mult = 0), limits = c(0, NA)) +
  scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
                  labels = scientific_10, name = "WT PTM\nabundance") +
  theme_classic(base_size = 6, base_family = "Arial") +
  theme(panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_line(linewidth = 0.2),
        legend.position = 'right', 
        legend.key = element_blank(),
        legend.key.size = unit(1.0, units = 'mm'),
        legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
        legend.direction = "vertical",
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.box.spacing = unit(-3, 'mm'),
        axis.text = element_text(colour = 'black'),
        axis.ticks.length = unit(1, units = 'mm'),
        axis.ticks = element_line(colour = 'black'),
        axis.title = element_text(size = 5),
        axis.line = element_line(linewidth = 0.2),
        plot.background = element_blank() ) -> p_KO

p_KO

Note: these WT and KO samples were processed at different moment, so there could be some batch effect difference that limits the statistical power.

Save to pdf.

Code
ggsave(path = pdf_dir_fig4, filename = "KO_vs_WT_Volcano_w_labels.pdf",
       plot = p_KO, device = cairo_pdf, units = "cm",
       height = 12, width = 10)

I could still add:

  • the rescues volcano plots vs the KO and KO+S vs KO+L

  • rescue efficiency

  • heatmap with all PTMs.

4 Supplementary Figure Panels

4.1 K27-K40 peptide modification ratios

4.1.1 H3.3K27

There’s no H3.3 signal in the main figure, here I plot only H3.3 signal.

Code
data_K27_36 |>
  group_by(Condition, K27, Histone) |>
  mutate(Average_Ratio = mean(Ratio)) |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3) |>
  subset(Histone %in% "H3.3") |>
  mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) ) |>
  ggplot() +
  aes(x = Condition, stratum = K27, alluvium = K27, 
      y = Ratio_Sum, fill = K27) +
  geom_alluvium(alpha = 0.5, width = 0.4) +
  geom_stratum2(alpha = 1, width = 0.6, colour = "black") +
  geom_fit_text(aes(label = K27 ),
                stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"),
                padding.y = grid::unit(0.2, "mm"),
                size = 7, family = "Arial", show.legend = F) +
  scale_fill_manual(values = c("K27ac" = "#feac81", "K27me1" ="#ced1af", 
                               "K27me2" = "#748f46", "K27me3" = "#47632a", 
                               "unmod" = "gray" ),
                    name = 'H3.3') +
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 
                              'KO+L', 'KO+S')) +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing 
  labs(y = "Average Ratio") +
  alluvial_thm -> p_Alluvial_H33K27

p_Alluvial_H33K27 

Save to pdf.

Code
ggsave(filename = "FigS4B_H3.3_only_K27_Mrgd_Alluvial.pdf", path = pdf_dir_fig4, 
       plot = p_Alluvial_H33K27, device = cairo_pdf, units = "cm",
       width = 6.5, height = 4.85)

4.1.2 H3 and H3.3 K36

Plot the pileup of H3 & H3.3 K36 methylation.

Code
data_K27_36 |>
  group_by(Condition, K36, Histone) |>
  mutate(Average_Ratio = mean(Ratio)) |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3) |>
  subset(Histone == c("H3", "H3.3") ) |>
  mutate(K36 = factor(K36, levels = c("K36me3", "K36me2", "K36me1", "unmod") ) ) |>
  ggplot() +
  aes(x = Condition, stratum = K36, alluvium = K36, 
      y = Ratio_Sum, fill = K36) +
  facet_wrap(~ Histone, scales = 'free_y') +
  geom_alluvium(alpha = 0.5, width = 0.4) +
  geom_stratum(alpha = 1, width = 0.6, linewidth = 0.1, colour = "black") +
  geom_fit_text(aes(label = K36 ), stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"), size = 7, family = "Arial",
                padding.y = grid::unit(0.2, "mm"),
                show.legend = F) +
  scale_fill_manual(values = c("K36me1" = "#87B7AA", "K36me2" = "#8FA9C2",
                               "K36me3" = "#8386B7", "unmod" = "gray" ) ) +
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 'KO+L', 'KO+S')
                   ) +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing 
  labs(y = "Average Ratio") +
  alluvial_thm +
  theme(legend.title = element_blank()) -> p_Alluvial_K36

p_Alluvial_K36

Save to pdf.

Code
ggsave(path = pdf_dir_fig4, filename = "FigS4C_H3_H3.3_K36_Mrgd_Alluvial.pdf",
       plot = p_Alluvial_K36, device = cairo_pdf, units = "cm",
       width = 12.0, height = 5)

4.2 Extra plots

Here I plot extra histone PTMs that were not included in the publication.

4.2.1 H3 K4

Plot H3K4 modifications.

Code
data_dia_auto |>
  group_by(Condition, Modification) |>
  subset(Peptide_Start >= 3 & Peptide_End <= 8) |>
  mutate(Average_Ratio = mean(Ratio))  |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
  ggplot() +
  aes(x = Condition, stratum = Modification, alluvium = Modification, 
      y = Ratio_Sum, fill = Modification) +
  geom_alluvium(alpha = 0.4, width = 0.4) +
  geom_stratum2(alpha = 1, width = 0.6, size = 0.1) +
  geom_fit_text(aes(label = Modification ),
                stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"),
                padding.y = grid::unit(0.2, "mm"),
                size = 7, family = "Arial", show.legend = F) +
  scale_fill_manual(values = met.brewer(name = "Demuth", direction = 1, n = 5), name = 'H3' ) +
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 
                              'KO+L', 'KO+S')) +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing 
  labs(y = "Average Ratio") +
  alluvial_thm +
  theme(legend.text = element_text(margin = margin(l = -9, unit = "mm")), 
        legend.key.width = unit(6.5, "mm") ) -> p_Alluvial_K4

p_Alluvial_K4

Save to pdf.

Code
ggsave(path = pdf_dir_fig4, filename = "H3K4_Mrgd_Alluvial.pdf",
       plot = p_Alluvial_K4, device = cairo_pdf, units = "cm",
       width = 6.5, height = 4.85)

4.2.2 H3 K9

Plot all combinatorial K9 and K14 modifications on histone H3.

Code
data_dia_auto |>
  group_by(Condition, Modification) |>
  subset(Histone == "H3") |>
  subset(Peptide_Start >= 9 & Peptide_End <= 17) |>
  mutate(Average_Ratio = mean(Ratio) ) |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
  ggplot() +
  aes(x = Condition, stratum = Modification, alluvium = Modification, 
      y = Ratio_Sum, fill = Modification) +
  geom_alluvium(alpha = 0.4, width = 0.4) +
  geom_stratum2(alpha = 1, width = 0.6, size = 0.1) +
  geom_fit_text(aes(label = Modification ),
                stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"),
                padding.y = grid::unit(0.2, "mm"),
                size = 7, family = "Arial", show.legend = F) +
  scale_fill_manual(values = met.brewer(name = "Hokusai1", direction = -1, n = 15), name = 'H3' ) + 
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 
                              'KO+L', 'KO+S')) +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  labs(y = "Average Ratio") +
  alluvial_thm +
  theme(legend.text = element_text(margin = margin(l = -12.0, unit = "mm")), 
        legend.key.width = unit(11.2, "mm"),
        legend.key.height = unit(1, 'mm'),
        legend.spacing.x = unit(1, 'mm') ) -> p_Alluvial_K9

p_Alluvial_K9

Code
ggsave(path = pdf_dir_fig4, filename = "H3K9_Mrgd_Alluvial.pdf",
       plot = p_Alluvial_K9, device = cairo_pdf, units = "cm",
       width = 6.4, height = 4.85)

4.2.3 All combinations of H3/H3.3 K27 K36

Plot all combinatorial K27 and K36 modifications on histone H3 and H3.3 variant.

Code
data_dia_auto |>
  group_by(Condition, Modification, Histone) |>
  subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
  mutate(Average_Ratio = mean(Ratio))  |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
  ggplot() +
  aes(x = Condition, stratum = Modification, alluvium = Modification, 
      y = Ratio_Sum, fill = Modification) +
  facet_wrap(~ Histone) +
  geom_alluvium(alpha = 0.4, width = 0.4) +
  geom_stratum(alpha = 1, width = 0.6, size = 0.1) +

  geom_fit_text(aes(label = Modification ),
                stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"),
                padding.y = grid::unit(0.2, "mm"),
                size = 7, family = "Arial", show.legend = F) +
  scale_fill_manual(values = met.brewer(name = "Hiroshige", direction = 1, n = 15) ) +
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 
                              'KO+L', 'KO+S'), name = 'H3 & H3.3') +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  labs(y = "Average Ratio") +
  alluvial_thm +
  theme(legend.text = element_text(margin = margin(l = -14.5, unit = "mm")), 
        legend.key.width = unit(13.75, "mm"),
        legend.key.height = unit(1, 'mm'),
        legend.spacing.x = unit(1, 'mm') ) -> p_Alluvial_K27_K36

p_Alluvial_K27_K36

This is the best representative picture of K27-K40 modifications relative abundance.

Code
ggsave(path = pdf_dir_fig4, filename = "H3_H3.3_K27_K36_Mrgd_Alluvial.pdf",
       plot = p_Alluvial_K27_K36, device = cairo_pdf, units = "cm",
       width = 14, height = 5)

End analysis.

5 Session Info

Code
sessioninfo::session_info()
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.1.2 (2021-11-01)
 os       macOS Catalina 10.15.7
 system   x86_64, darwin17.0
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Europe/Madrid
 date     2023-09-04
 pandoc   3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 ! package              * version   date (UTC) lib source
   ade4                   1.7-22    2023-02-06 [1] CRAN (R 4.1.2)
   affy                   1.72.0    2021-10-26 [1] Bioconductor
   affyio                 1.64.0    2021-10-26 [1] Bioconductor
   annotate               1.72.0    2021-10-26 [1] Bioconductor
   AnnotationDbi          1.56.2    2021-11-09 [1] Bioconductor
   assertthat             0.2.1     2019-03-21 [1] CRAN (R 4.1.0)
   Biobase              * 2.54.0    2021-10-26 [1] Bioconductor
   BiocFileCache          2.2.1     2022-01-23 [1] Bioconductor
   BiocGenerics         * 0.40.0    2021-10-26 [1] Bioconductor
   BiocManager            1.30.22   2023-08-08 [1] CRAN (R 4.1.2)
   BiocParallel           1.28.3    2021-12-09 [1] Bioconductor
   biomaRt                2.50.3    2022-02-06 [1] Bioconductor
   Biostrings             2.62.0    2021-10-26 [1] Bioconductor
   bit                    4.0.5     2022-11-15 [1] CRAN (R 4.1.2)
   bit64                  4.0.5     2020-08-30 [1] CRAN (R 4.1.0)
   bitops                 1.0-7     2021-04-24 [1] CRAN (R 4.1.0)
   blob                   1.2.4     2023-03-17 [1] CRAN (R 4.1.2)
   bslib                  0.5.1     2023-08-11 [1] CRAN (R 4.1.2)
   cachem                 1.0.8     2023-05-01 [1] CRAN (R 4.1.2)
   Cairo                * 1.6-1     2023-08-18 [1] CRAN (R 4.1.2)
   callr                  3.7.3     2022-11-02 [1] CRAN (R 4.1.2)
   cellranger             1.1.0     2016-07-27 [1] CRAN (R 4.1.0)
   circlize               0.4.15    2022-05-10 [1] CRAN (R 4.1.2)
   cli                    3.6.1     2023-03-23 [1] CRAN (R 4.1.2)
   clue                   0.3-64    2023-01-31 [1] CRAN (R 4.1.2)
   cluster                2.1.4     2022-08-22 [1] CRAN (R 4.1.2)
   codetools              0.2-19    2023-02-01 [1] CRAN (R 4.1.2)
   colorspace             2.1-0     2023-01-23 [1] CRAN (R 4.1.2)
   ComplexHeatmap         2.10.0    2021-10-26 [1] Bioconductor
   crayon                 1.5.2     2022-09-29 [1] CRAN (R 4.1.2)
   crosstalk              1.2.0     2021-11-04 [1] CRAN (R 4.1.0)
   csaw                 * 1.28.0    2021-10-26 [1] Bioconductor
   curl                   5.0.2     2023-08-14 [1] CRAN (R 4.1.2)
   DBI                    1.1.3     2022-06-18 [1] CRAN (R 4.1.2)
   dbplyr                 2.3.3     2023-07-07 [1] CRAN (R 4.1.2)
   DelayedArray           0.20.0    2021-10-26 [1] Bioconductor
   DEP                  * 1.16.0    2021-10-26 [1] Bioconductor
   desc                   1.4.2     2022-09-08 [1] CRAN (R 4.1.2)
   DESeq2               * 1.34.0    2021-10-26 [1] Bioconductor
   devtools               2.4.5     2022-10-11 [1] CRAN (R 4.1.2)
   digest                 0.6.33    2023-07-07 [1] CRAN (R 4.1.2)
   doParallel             1.0.17    2022-02-07 [1] CRAN (R 4.1.2)
   dplyr                * 1.1.2     2023-04-20 [1] CRAN (R 4.1.2)
   DT                   * 0.29      2023-08-29 [1] CRAN (R 4.1.2)
   edgeR                  3.36.0    2021-10-26 [1] Bioconductor
   ellipsis               0.3.2     2021-04-29 [1] CRAN (R 4.1.0)
   evaluate               0.21      2023-05-05 [1] CRAN (R 4.1.2)
   fansi                  1.0.4     2023-01-22 [1] CRAN (R 4.1.2)
   farver                 2.1.1     2022-07-06 [1] CRAN (R 4.1.2)
   fastmap                1.1.1     2023-02-24 [1] CRAN (R 4.1.2)
   fdrtool                1.2.17    2021-11-13 [1] CRAN (R 4.1.0)
   filelock               1.0.2     2018-10-05 [1] CRAN (R 4.1.0)
   forcats                1.0.0     2023-01-29 [1] CRAN (R 4.1.2)
   foreach                1.5.2     2022-02-02 [1] CRAN (R 4.1.2)
   foreign                0.8-84    2022-12-06 [1] CRAN (R 4.1.2)
   fs                     1.6.3     2023-07-20 [1] CRAN (R 4.1.2)
   genefilter             1.76.0    2021-10-26 [1] Bioconductor
   geneplotter            1.72.0    2021-10-26 [1] Bioconductor
   generics               0.1.3     2022-07-05 [1] CRAN (R 4.1.2)
   GenomeInfoDb         * 1.30.1    2022-01-30 [1] Bioconductor
   GenomeInfoDbData       1.2.7     2023-08-29 [1] Bioconductor
   GenomicRanges        * 1.46.1    2021-11-18 [1] Bioconductor
   GetoptLong             1.0.5     2020-12-15 [1] CRAN (R 4.1.0)
   ggalluvial             0.12.5    2023-02-22 [1] CRAN (R 4.1.2)
   ggfittext              0.10.0    2023-04-04 [1] CRAN (R 4.1.2)
   ggplot2                3.4.3     2023-08-14 [1] CRAN (R 4.1.2)
   ggrepel                0.9.3     2023-02-03 [1] CRAN (R 4.1.2)
   ggseqlogo              0.1       2017-07-25 [1] CRAN (R 4.1.0)
   GlobalOptions          0.1.2     2020-06-10 [1] CRAN (R 4.1.0)
   glue                   1.6.2     2022-02-24 [1] CRAN (R 4.1.2)
   gmm                    1.8       2023-06-06 [1] CRAN (R 4.1.2)
   gtable                 0.3.4     2023-08-21 [1] CRAN (R 4.1.2)
   hexbin                 1.28.3    2023-03-21 [1] CRAN (R 4.1.2)
   hms                    1.1.3     2023-03-21 [1] CRAN (R 4.1.2)
   htmltools              0.5.6     2023-08-10 [1] CRAN (R 4.1.2)
   htmlwidgets            1.6.2     2023-03-17 [1] CRAN (R 4.1.2)
   httpuv                 1.6.11    2023-05-11 [1] CRAN (R 4.1.2)
   httr                   1.4.7     2023-08-15 [1] CRAN (R 4.1.2)
   impute                 1.68.0    2021-10-26 [1] Bioconductor
   imputeLCMD             2.1       2022-06-10 [1] CRAN (R 4.1.2)
   IRanges              * 2.28.0    2021-10-26 [1] Bioconductor
   iterators              1.0.14    2022-02-05 [1] CRAN (R 4.1.2)
   jquerylib              0.1.4     2021-04-26 [1] CRAN (R 4.1.0)
   jsonlite               1.8.7     2023-06-29 [1] CRAN (R 4.1.2)
   KEGGREST               1.34.0    2021-10-26 [1] Bioconductor
   knitr                  1.43      2023-05-25 [1] CRAN (R 4.1.2)
   labeling               0.4.2     2020-10-20 [1] CRAN (R 4.1.0)
   later                  1.3.1     2023-05-02 [1] CRAN (R 4.1.2)
   lattice                0.21-8    2023-04-05 [1] CRAN (R 4.1.2)
   lifecycle              1.0.3     2022-10-07 [1] CRAN (R 4.1.2)
   limma                  3.50.3    2022-04-07 [1] Bioconductor
   locfit                 1.5-9.8   2023-06-11 [1] CRAN (R 4.1.2)
   magick                 2.7.5     2023-08-07 [1] CRAN (R 4.1.2)
   magrittr               2.0.3     2022-03-30 [1] CRAN (R 4.1.2)
   MALDIquant             1.22.1    2023-03-20 [1] CRAN (R 4.1.2)
   MASS                   7.3-60    2023-05-04 [1] CRAN (R 4.1.2)
   Matrix                 1.6-1     2023-08-14 [1] CRAN (R 4.1.2)
   MatrixGenerics       * 1.6.0     2021-10-26 [1] Bioconductor
   matrixStats          * 1.0.0     2023-06-02 [1] CRAN (R 4.1.2)
   memoise                2.0.1     2021-11-26 [1] CRAN (R 4.1.0)
   metapod                1.2.0     2021-10-26 [1] Bioconductor
   MetBrewer              0.2.0     2022-03-21 [1] CRAN (R 4.1.2)
   mime                   0.12      2021-09-28 [1] CRAN (R 4.1.0)
   miniUI                 0.1.1.1   2018-05-18 [1] CRAN (R 4.1.0)
   msa                    1.26.0    2021-10-26 [1] Bioconductor
   MsCoreUtils            1.6.2     2022-02-24 [1] Bioconductor
   MSnbase                2.20.4    2022-01-16 [1] Bioconductor
   munsell                0.5.0     2018-06-12 [1] CRAN (R 4.1.0)
   mvtnorm                1.2-3     2023-08-25 [1] CRAN (R 4.1.2)
   mzID                   1.32.0    2021-10-26 [1] Bioconductor
   mzR                    2.28.0    2021-10-27 [1] Bioconductor
   ncdf4                  1.21      2023-01-07 [1] CRAN (R 4.1.2)
 R niar                 * 0.3.0     <NA>       [?] <NA>
   norm                   1.0-11.1  2023-06-18 [1] CRAN (R 4.1.2)
   patchwork              1.1.3     2023-08-14 [1] CRAN (R 4.1.2)
   pcaMethods             1.86.0    2021-10-26 [1] Bioconductor
   pheatmap             * 1.0.12    2019-01-04 [1] CRAN (R 4.1.0)
   pillar                 1.9.0     2023-03-22 [1] CRAN (R 4.1.2)
   pkgbuild               1.4.2     2023-06-26 [1] CRAN (R 4.1.2)
   pkgconfig              2.0.3     2019-09-22 [1] CRAN (R 4.1.0)
   pkgload                1.3.2.1   2023-07-08 [1] CRAN (R 4.1.2)
   plyr                   1.8.8     2022-11-11 [1] CRAN (R 4.1.2)
   png                    0.1-8     2022-11-29 [1] CRAN (R 4.1.2)
   preprocessCore         1.56.0    2021-10-26 [1] Bioconductor
   prettyunits            1.1.1     2020-01-24 [1] CRAN (R 4.1.0)
   processx               3.8.2     2023-06-30 [1] CRAN (R 4.1.2)
   profvis                0.3.8     2023-05-02 [1] CRAN (R 4.1.2)
   progress               1.2.2     2019-05-16 [1] CRAN (R 4.1.0)
   promises               1.2.1     2023-08-10 [1] CRAN (R 4.1.2)
   ProtGenerics           1.26.0    2021-10-26 [1] Bioconductor
   ps                     1.7.5     2023-04-18 [1] CRAN (R 4.1.2)
   psychTools             2.3.6     2023-06-18 [1] CRAN (R 4.1.2)
   purrr                  1.0.2     2023-08-10 [1] CRAN (R 4.1.2)
   R6                     2.5.1     2021-08-19 [1] CRAN (R 4.1.0)
   rappdirs               0.3.3     2021-01-31 [1] CRAN (R 4.1.0)
   RColorBrewer           1.1-3     2022-04-03 [1] CRAN (R 4.1.2)
   Rcpp                   1.0.11    2023-07-06 [1] CRAN (R 4.1.2)
   RCurl                  1.98-1.12 2023-03-27 [1] CRAN (R 4.1.2)
   readr                * 2.1.4     2023-02-10 [1] CRAN (R 4.1.2)
   readxl               * 1.4.3     2023-07-06 [1] CRAN (R 4.1.2)
   remotes                2.4.2.1   2023-07-18 [1] CRAN (R 4.1.2)
   rjson                  0.2.21    2022-01-09 [1] CRAN (R 4.1.2)
   rlang                  1.1.1     2023-04-28 [1] CRAN (R 4.1.2)
   rmarkdown              2.24      2023-08-14 [1] CRAN (R 4.1.2)
   rprojroot              2.0.3     2022-04-02 [1] CRAN (R 4.1.2)
   Rsamtools              2.10.0    2021-10-26 [1] Bioconductor
   RSQLite                2.3.1     2023-04-03 [1] CRAN (R 4.1.2)
   rstudioapi             0.15.0    2023-07-07 [1] CRAN (R 4.1.2)
   S4Vectors            * 0.32.4    2022-03-29 [1] Bioconductor
   sandwich               3.0-2     2022-06-15 [1] CRAN (R 4.1.2)
   sass                   0.4.7     2023-07-15 [1] CRAN (R 4.1.2)
   scales                 1.2.1     2022-08-20 [1] CRAN (R 4.1.2)
   seqinr                 4.2-30    2023-04-05 [1] CRAN (R 4.1.2)
   sessioninfo            1.2.2     2021-12-06 [1] CRAN (R 4.1.0)
   shape                  1.4.6     2021-05-19 [1] CRAN (R 4.1.0)
   shiny                  1.7.5     2023-08-12 [1] CRAN (R 4.1.2)
   shinydashboard         0.7.2     2021-09-30 [1] CRAN (R 4.1.0)
   stringi                1.7.12    2023-01-11 [1] CRAN (R 4.1.2)
   stringr                1.5.0     2022-12-02 [1] CRAN (R 4.1.2)
   SummarizedExperiment * 1.24.0    2021-10-26 [1] Bioconductor
   survival               3.5-7     2023-08-14 [1] CRAN (R 4.1.2)
   tibble                 3.2.1     2023-03-20 [1] CRAN (R 4.1.2)
   tidyr                  1.3.0     2023-01-24 [1] CRAN (R 4.1.2)
   tidyselect             1.2.0     2022-10-10 [1] CRAN (R 4.1.2)
   tmvtnorm               1.5       2022-03-22 [1] CRAN (R 4.1.2)
   tzdb                   0.4.0     2023-05-12 [1] CRAN (R 4.1.2)
   urlchecker             1.0.1     2021-11-30 [1] CRAN (R 4.1.0)
   usethis                2.2.2     2023-07-06 [1] CRAN (R 4.1.2)
   utf8                   1.2.3     2023-01-31 [1] CRAN (R 4.1.2)
   vctrs                  0.6.3     2023-06-14 [1] CRAN (R 4.1.2)
   vroom                  1.6.3     2023-04-28 [1] CRAN (R 4.1.2)
   vsn                    3.62.0    2021-10-26 [1] Bioconductor
   withr                  2.5.0     2022-03-03 [1] CRAN (R 4.1.2)
   xfun                   0.40      2023-08-09 [1] CRAN (R 4.1.2)
   XICOR                  0.4.1     2023-04-21 [1] CRAN (R 4.1.2)
   XML                    3.99-0.14 2023-03-19 [1] CRAN (R 4.1.2)
   xml2                   1.3.5     2023-07-06 [1] CRAN (R 4.1.2)
   xtable                 1.8-4     2019-04-21 [1] CRAN (R 4.1.0)
   XVector                0.34.0    2021-10-26 [1] Bioconductor
   yaml                   2.3.7     2023-01-23 [1] CRAN (R 4.1.2)
   zlibbioc               1.40.0    2021-10-26 [1] Bioconductor
   zoo                    1.8-12    2023-04-13 [1] CRAN (R 4.1.2)

 [1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library

 R ── Package was removed from disk.

──────────────────────────────────────────────────────────────────────────────

References

Yuan, Zuo-Fei, Simone Sidoli, Dylan M Marchione, Johayra Simithy, Kevin A Janssen, Mary R Szurgot, and Benjamin A Garcia. 2018. EpiProfile 2.0: A Computational Platform for Processing Epi-Proteomics Mass Spectrometry Data.” Journal of Proteome Research 17 (7): 2533–41. https://doi.org/10.1021/acs.jproteome.8b00133.